/* Part of the code is copied form test_s2_semi_memo_for.c */
#include <errno.h>
#include <cmath>
#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <algorithm>
#include <string>

#include "fftw3.h"
#include "makeweights.h"
#include "cospmls.h"
#include "FST_semi_memo.h"
#include "des_gen.h"

using namespace std;

vector< vector<double> > SphHarmo :: normalize(vector< vector<double> > undes) {
	vector<vector<double> > result;
	vector<double> min,max;
	vector<double> temp;
	double temp_value;
	double min_value, max_value;
	
    for (int i = 0; i < undes.size(); i++) {
		min.push_back(*min_element(undes[i].begin(),undes[i].end()));
		max.push_back(*max_element(undes[i].begin(),undes[i].end()));
	}

	min_value = *min_element(min.begin(),min.end());
	max_value = *max_element(max.begin(),max.end());
	
    for(int i = 0; i < undes.size(); i++) {
		temp.clear();
		for(int j = 0; j < undes[i].size(); j++) {
			temp_value = (undes[i][j]-min_value) / (max_value-min_value) * 10;
			temp.push_back(temp_value);
		}
		result.push_back(temp);
	}
	return result;
}

void SphHarmo :: getdes(const char* filename,int BW) {
	int size;
	int dummy, rtemp;
	int cutoff, order;
	int rank, howmany_rank;
	vector< vector<double> > des; // 32*BW
	vector< vector<double> > undes(32,vector<double>(BW));
	double *rdata;
	double *rcoeffs, *icoeffs;
	double *weights;
	double *seminaive_naive_tablespace, *workspace;
	double **seminaive_naive_table;
	fftw_plan dctPlan, fftPlan;
	fftw_iodim dims[1], howmany_dims[1];
	string buffer;

	cutoff = BW;
	size = 2*BW;
	double idata[(size*size)];
	fill(idata,idata+size*size,0.0);

	/* allocate memory */
	rdata = (double *) malloc(sizeof(double) * (size * size));
	rcoeffs = (double *) malloc(sizeof(double) * (BW * BW));
	icoeffs = (double *) malloc(sizeof(double) * (BW * BW));
	weights = (double *) malloc(sizeof(double) * 4 * BW);
	seminaive_naive_tablespace =
				(double *) malloc(sizeof(double) *
				(Reduced_Naive_TableSize(BW,cutoff) +
				Reduced_SpharmonicTableSize(BW,cutoff)));
	workspace = (double *) malloc(sizeof(double)*((8 * (BW*BW))+(7 * BW)));

	/* check to see if all the memory has been allocated */
	if ( (rdata == NULL) ||
       (rcoeffs == NULL) || (icoeffs == NULL) ||
       (seminaive_naive_tablespace == NULL) ||
       (workspace == NULL) ) {
		perror("Error in allocating memory");
		exit(1);
	}

	/* now precompute the Legendres */
	/* Generating seminaive_naive tables */
	seminaive_naive_table = SemiNaive_Naive_Pml_Table(BW, cutoff,
						    seminaive_naive_tablespace,
						    workspace);

	/* construct fftw plans */
	/* make DCT plan -> using the GURU interface to
	 * execute these plans within the routines*/
	/* forward DCT */
	dctPlan = fftw_plan_r2r_1d( 2*BW, weights, rdata,
			      FFTW_REDFT10, FFTW_ESTIMATE );
      
	/* fftw "preamble";
     note that this plan places the output in a transposed array */
	rank = 1;
	dims[0].n = 2*BW;
	dims[0].is = 1;
	dims[0].os = 2*BW;
	howmany_rank = 1;
	howmany_dims[0].n = 2*BW;
	howmany_dims[0].is = 2*BW;
	howmany_dims[0].os = 1;
  
	/* forward fft */
	fftPlan = fftw_plan_guru_split_dft( rank, dims,
				      howmany_rank, howmany_dims,
				      rdata, idata,
				      workspace, workspace+(4*BW*BW),
				      FFTW_ESTIMATE );
	/* now make the weights */
	makeweights( BW, weights );
	
	string temp_fn(filename);
	string::size_type nbegin;
	string::size_type nend;

	nbegin = temp_fn.rfind('/',temp_fn.size());
	nend = temp_fn.rfind('.',temp_fn.size());
	if (nbegin == string::npos || nend == string::npos) {
		cout<<"wrong file name"<<endl;
		exit(0);
	}
	string objname = temp_fn.substr((nbegin+1),(nend-nbegin-1));
	string dataname = "./data/"+objname+".dat";
	ifstream ifp (dataname.c_str());

	for (int rindex = 0; rindex < 32; rindex++) {
		getline(ifp,buffer);
		for (int i = 0; i < size*size; i++){
			getline(ifp,buffer);
			stringstream ss(buffer);
			ss >> rdata[i];
		}
		/* now do the forward spherical transform */
		FST_semi_memo(rdata, idata, rcoeffs, icoeffs, BW,
				seminaive_naive_table, workspace,0 ,cutoff,
				&dctPlan, &fftPlan, weights );

		/* now write out coefficients*/
		double temp = 0;
		/* code format */
		for (int l = 0; l < BW; l++) {
			temp = 0;
			for (int m = -l; m < l+1; m++) {
				dummy = seanindex(m, l, BW);
				temp += (pow(rcoeffs[dummy],2)+pow(icoeffs[dummy],2));
			}
			temp = sqrt(temp);
			undes[rindex][l] = temp;
		}
	}

	ifp.close();
	des = normalize(undes);
	
	string ofilename = "out/" + objname + ".dat";
	ofstream ofp(ofilename.c_str());
	for (int i = 1; i < 33; i++) {
		ofp << "# " << i << endl;
		for (int j = 0; j < BW; j++) {
			ofp << des[i-1][j] << endl;
		}
	}
	ofp.close();

	/* clean up */
	fftw_destroy_plan(fftPlan);
	fftw_destroy_plan(dctPlan);
	free(workspace);
	free(seminaive_naive_table);
	free(seminaive_naive_tablespace);
	free(weights);
	free(icoeffs);
	free(rcoeffs);
	free(rdata);

	cout << "# " << objname << ".dat " << "is added" << endl;
	cout << "------------------------------------------" << endl;
}
